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The instability of r-mode oscillations in rapidly rotating neutron stars has attracted attention as a potential 
mechanism for producing high frequency, almost periodic gravitational waves. The analyses carried so far have 
shown the existence of these modes and have considered damping by shear and bulk viscosity. However, the 
magnetohydrodynamic coupling of the modes with a stellar magnetic field and its role in the damping of the 
instability has not been fully investigated yet. Following our introductory paper flip, we here discuss in more 
detail the existence of secular higher-order kinematical effects which will produce toroidal fluid drifts. We 
also define the sets of equations that account for the time evolution of the magnetic fields produced by these 
secular velocity fields and show that the magnetic fields produced can reach equipartition in less than a year. 
The full numerical calculations as well as the evaluation of the impact of strong magnetic fields on the onset and 
evolution of the r-mode instability will be presented in a companion paper Pp. 
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|^ ■ I. INTRODUCTION 

o. 

The basic properties of r-mode oscillations in Newtonian rotating stars were investigated in 1978 by Papaloizou and Pringle [|[j 
in an attempt to explain the short-period oscillations seen in cataclysmic variables in terms of non-radial oscillations of rotating 
white dwarfs. Subsequent work [£Q] has increased our understanding of these oscillations, which have close similarities with the 
Rossby waves observed in the Earth's atmosphere and oceans. Almost twenty years later, r-mode oscillations have become the 
focus of renewed attention when they were shown to be unstable to the emission of gravitational radiation. The first calculations 
in this sense were carried out by Andersson [g] and by Friedman and Morsink [[7]]. Since then, the literature on the subject 
>• has been growing rapidly. Exhaustive reviews of our present understanding and of unresolved issues concerning the r-mode 
J^j ' instability can be found in [g|] and in [g] . 

A significant difference from previously investigated mode-instabilities is that gravitational radiation couples with r-mode 
C^ ■ oscillations primarily through time-varying mass-current multipole moments rather than through the usual time-varying mass 
multipole moments. Coherent, large-scale fluid currents in a hot plasma with high electrical conductivity, such as neutron star 
matter, represent the basic conditions for the generation of large scale magnetic fields. This paper is devoted to studying whether 
intense magnetic fields can be produced as a result of the onset and growth of the r-mode instability. In particular, we here 
extend the discussion first presented in ||l|] about the existence of secular velocity fields which interact with the magnetic field 
pre-existing in the neutron star. We show in more detail that such secular effects can be derived with confidence from linearized 
expressions and are responsible for differential rotation, both in the radial and in the polar directions. As a result, secular toroidal 
drifts appear on isobaric surfaces and generate large-scale magnetic fields. We also derive here the sets of equations necessary 
for the calculation of the magnetic fields produced in this way. In a companion paper [Efl (hereafter paper II) we will present in 
detail the results of the numerical solution of the equations presented here and comment on the impact that the large magnetic 
fields produced will have on the onset and development of the instability. 

The organization of the paper is as follows: Section y introduces the physical framework of our approach and synthesizes the 
results. Section pH, takes a closer look at the linear equations of motion for fiducial fluid elements on isobaric surfaces. Starting 
from those we then derive the higher-order secular velocity field and comment on the validity of our results. Section |v] is devoted 
to the formulation of the set of equations necessary for the calculation of the magnetic field produced by the secular velocity 
field derived in the previous Section. In particular, we discuss two different approaches and comment on the corresponding 
advantages and shortcomings. Section M finally presents our conclusions and refers the reader to the numerical results and their 
astrophysical consequences which are presented in paper II. 



II. PHYSICAL PICTURE 

The energy budget governing the evolution of the r-mode instability is traditionally assumed to be regulated only by the emis- 
sion of gravitational waves and viscosity, which act as sources and sinks of energy, respectively [I1Q-IT2] (In a frame corotating 
with the star, the emission of gravitational waves is seen as an input of energy for the mode.)- However, it is natural to expect that 
other sources and sinks should intervene in this budget, most notably the loss of rotational and mode energy to electromagnetic 
radiation [ |l3| ] and to the coupling of the mass-currents produced by the oscillations with the magnetic field present in the neutron 
star [jj], In the case of the r-mode instability, this is particularly relevant. The reason for this is twofold: firstly, because the 
oscillations produce large scale mass-currents and, secondly, because the generation of magnetic field is a generic feature of 
shearing flows perpendicular to magnetic field lines in a highly conducting plasma, such as hot neutron star matter. This reflects 
the fact that, in such conditions, the magnetic field is predominantly advected with the fluid. When the electrical conductivity is 
infinite (the ideal MHD limit), magnetic field lines are "frozen" in the fluid and move entirely with it (Alfven Theorem, fll5|]). As 
pointed out by Spruit [Hq], the generation of magnetic fields under these circumstances can be so efficient that extremely intense 
magnetic fields could be created during the instability. When these become buoyant-unstable they can then generate powerful 
flashes of 7-rays. 

This paper, similarly to the one preceding it [[!]] and the one following [Q], aims at determining the strength of this coupling 
and the consequences that it will have on the onset and evolution of the instability. In particular, we will show that the kinematic 
properties of the r-mode oscillations will give rise to a secular velocity field which, once coupled to any seed magnetic field, 
will produce exponentially growing magnetic fields as the instability develops. A detailed analysis of this coupling, which 
inevitably introduces a change in the character of the modes, is beyond the scope of the present paper. However, the important 
qualitative effects can be estimated rather simply. If the magnetic field produced as a result of this coupling is strong enough, 
it will significantly distort the r-mode oscillations so as prevent the amplification of r modes by gravitational radiation. If the 
magnetic field is initially weak, it will be subsequently amplified and cause the instability to die out as the star spins down. 

III. KINEMATIC PROPERTIES OF THE LINEAR R MODES 

As mentioned above, a key role in the coupling between r-mode oscillations and the magnetic field is played by the kinematic 
properties of the modes. As we will discuss below, a careful look at the equations of motion will reveal nonlinear effects which 
manifest themselves mostly in secular velocity fields. Within this Section we will assume that the r-mode oscillations have 



amplitudes which are constant in time, that the nonlinear coupling among different modes is negligible [ 17 1, and that there is no 
magnetic field. The analysis of the kinematical properties of r-mode oscillations in the case of a time-growing mode amplitude 
will be discussed in paper II (cf. Section II. A) 

A. The Eulerian perturbation velocity field 

We begin our analysis by considering the motion of fiducial fluid elements on an isobaric surface of a rotating star experiencing 
r-mode oscillations. For a Newtonian inviscid star, these are solutions of the perturbed hydrodynamic equations having Eulerian 
velocity perturbations of "axial type" |]10|,[7||. In an orthonormal basis, and at first order in the star's unperturbed angular velocity 
f2, such perturbations may be written as 

&v 1 (r,e,ht) = otoR(^ t Yf m e i ° t , (1) 

where R is the radius of the unperturbed star and a is the frequency of the mode in the inertial frame. The dimensionless 
coefficient a describes the amplitude of the perturbation and we have indicated with and index "1" the velocity perturbations 
which are linear in a. In Eq. (fil), Yf m is the magnetic-type vector spherical harmonic, and may be defined in terms of the 
spherical harmonic functions Yg m by (see, e.g., [|l8||) 

[rV x (r\7Y (m )] . (2) 
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Consider now a frame instantaneously corotating with the star. In this frame, the differential equations governing the motion 
of fiducial fluid elements of an £ = m mode[]in the coordinate basis (i, r, 9, cf>) are 



1 Hereafter we will always refer to modes for which m — I. 
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and the dot refers to the total derivative with respect to the time coordinate. The angular frequency a of an r-mode in the inertial 
frame can be related to its_ angular frequency u in the corotating frame and the stellar angular velocity. At the lowest order in 9, 
this relation is given by | 



a = u) — £9 = 
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Throughout this paper and in the companion paper II we will restrict ourselves to the study of the r-mode instability in the 
slow rotation approximation, retaining only the the lowest order term in 9. Nevertheless, we will often consider neutron stars 
spinning very near the break-up limit. While this approximation might be a reasonable one [Il9lj, it is important to bear in mind 
that significant modifications could appear in our picture of the r-mode instability if the general relativistic rotational effects are 
fully taken into account. 



B. Nonlinear Motions of Fluid Elements with constant amplitude 



As will become apparent in Section IV , the linear-order equations (Bj) do not generate a significant secular magnetic field 
since they lead to unitary strain tensors [cf. Eq. (Jl4|)]. However, equations (||) can provide important information about the 
nonlinear motions of fluid elements and, in particular, about whether they lead to a secular drift velocity. When the nonlinear 
expressions are not available, in fact, a rather standard technique pC| , pl] ] allows to calculate second-order quantities from linear 
results. This is an approximation but in some relevant examples, such as sound waves and shallow water waves, one finds that 
the lowest-order nonlinear corrections to the linear velocity field make no contribution to the estimated velocity field: i.e. the 
drift velocity is given exactly to 0(a 2 ) by the linear velocity field. 

We have made use of this technique and obtained analytical expressions for the values of the velocity perturbations at second- 
order in a. In particular, we have expanded the equations of motion in powers of a, averaged over a gyration, and retained only 
the lowest-order non-vanishing term (see Appendix W for details). Interestingly, we find that a second-order in a drift velocity 
exists in the 0— direction and the total displacement in <fi from the onset of the oscillation at to up to time t is then found to be 
(see Appendix [ 
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with i 3> At. The matching coefficient np{6) is introduced to relate the instantaneous and secular velocities and is dependent on 
the mode number £ and on the 6 position on an isobaric surface of the star. For the I = 2, K2(0) = (l/2) 7 (5!/7r)(sin 2 6—2 cos 2 9) 
and hence the net displacement in the azimuthal direction after an oscillation is approximately 2na times the radius of gyration. 
The presence of a secular motion in Eq. (g|) is indicated by a non-periodic argument for the integral in the second expression on 
the left-hand-side. Note also that, for a star with constant angular velocity, a net secular azimuthal motion is obtained both when 
the mode's amplitude is constant [in which case this is 0(a 2 t)] and when the mode's amplitude is exponentially growing in time 
[in which case this is ©(a 2 , exp(2t/r GW )), where ao is the initial mode amplitude and r GW is the timescale for the instability to 
develop]. As we will discuss in paper II, this latter case will produce an even stronger azimuthal drift. 



At second-order in a there is no secular motion in the #— direction. 



The azimuthal drift velocity of a given fluid element is readily calculated from (g|). In an orthonormal basis and for the £ 
mode, this is 



v d (r,«,()^ K2 («)« 2 (f)fi(t)^)' 



e; 



(7) 



where ei is the unit vector in the cf> direction. It is worth underlining that the secular velocity field (Uh is responsible for a 
differential rotation in both the radial and polar directions. 

It is important to emphasize that using equations (Bj) to compute the displacement of an element of fluid expanding 9 and <j) 
in powers of a is not equivalent to considering nonlinear effects in the fluid equations. The order in a to which we compute 
the fluid displacement given by the linear velocity field (Bh is therefore a distinct (although related) question. We can compute 
the fluid motions predicted by the linear velocity field to any order we like, but if we compute them to orders higher than a 
it does not mean we have properly considered nonlinear fluid effects. However, using as guides analogous fluid-dynamical 
processes whose nonlinear behaviour is known, we expect the existence of a secular drift velocity of 0(a 2 ). Moreover, we 
expect the drift of fluid elements given by the velocity field v<j to be qualitatively correct and perhaps exact to 0(a 2 ). After 
this prediction was first mad e [U n, the existence of an 0(a 2 ) drift velocity and differential rotation has been verified both in 
analytical simplified models [g2|] and, at least qualitatively, through nonlinear numerical simulations p3|^4| ]. It should also be 
noted that the differential rotation described by Eq. O) is of kinematic nature and is set up on the timescale given by the rotation 
period of the star. This is much shorter than the timescale of gravitational wave emission over which the differential rotation, 
first proposed by Spruit [n60, is produced. 
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FIG. 1. Projected trajectories 8(t) sin 6(t) cos 4>(i), and <j>(t) sin 8(t) cos (j>(t) of seven fiducial fluid elements as seen in the corotating 
frame. The left diagram shows the projected trajectories over one period of fiducial fluid elements subject to an r-mode oscillation with £ = 2, 
a = 0.1, and SI = il B . All of the fluid elements are initially positioned on the 4>o = meridian and have different latitudes (indicated with 
stars). Continuous lines show the simultaneous positions of the fluid elements and different lines refer to different times. The trajectories of 
some representative fluid elements are shown with dashed circles. The small inset in the left diagram of Fig. Ill shows the initial positions 
(dashed line) and final positions (continuous line) of the different fluid elements. The right panel uses the same quantities as in the left one, 
but shows the trajectories of eight fluid elements for five oscillations. 



In Fig. fi] we show numerical integrations of the equations of motion (Bbh and ( pq ) on the northern hemisphere of the rotating 
star. In particular, the left panel of Fig. [I] shows the projected (9, <fi) trajectories over one period for fiducial fluid elements subject 
to an r-mode oscillation with £ = 2 (the coordinates are those of a reference frame corotating with the star). All of the fluid 
elements are initially positioned on the <j>o = meridian, but at different latitudes and these positions are indicated with stars 
(because of the polar dependence of equations [5H| and [[3c]], it is sufficient to study them in the northern hemisphere only). The 
continuous lines show the simultaneous positions of fluid elements at different times (some of these positions are indicated with 



circles), and the dashed lines are used to trace the trajectories of representative fluid elements. In solving equations ( pq) an d ( pq ) 
we have assumed that a has the constant value 0.1 and that the star is rotating at the break-up limit B = (2/3)^/irGp, with 
p = 3M/(AttR 3 ) being the average mass density. 

Several interesting properties can be deduced from Fig. |} a) the trajectories over one period uj/2tt are not closed; b) while the 
excursions in the 8 and <p directions are comparable, the net displacement Ac/) in longitude after one cycle is much larger than the 
net displacement A8 in latitude (indeed one can show that this is zero at 0(a 2 )); c) both Acf> and A9 are linearly proportional 
to a; d) the motions have a simple dependence on latitude. The small inset in the left diagram of Fig. [I] shows schematically the 
position of the fluid elements at the beginning of the oscillation (dashed line) and after one cycle (continuous line). The right 
panel of Fig. n] shows the trajectories of fiducial fluid elements during five full oscillations and the net displacement produced by 
their motions. 

Despite the complexity of the different trajectories, all of the information about their secular behaviour is contained in (H) and 
(Q). These expressions are at the basis of our study of the interaction between r-mode oscillations and any magnetic field that is 
initially present in the rotating neutron star. 

IV. EVOLUTION EQUATIONS FOR THE MAGNETIC FIELD 



In this Section we first discuss the model assumed for a magnetic neutron star (Section IV A ) and then illustrate the techniques 
employed to compute the evolution of an arbitrarily small seed magnetic field as a result of the secular drift produced by the 
r-mode oscillations. In particular, we have considered both Lagrangian and Eulerian formulations of the induction equation. 
The first exploits the possibility of calculating the strain tensor produced by the oscillation. The second approach, on the other 
hand, is more commonly used and treats the induction equation as a set of partial differential equations. For compactness, the 
numerical results obtained using both of these methods are presented in Paper II. 

A. Simplified model of magnetic neutron stars 

The properties of the magnetic fields of neutron stars still contain a number of unresolved aspects, mostly connected with the 



properties of matter at very high densities [25 ]. For simplicity, we will assume that the stellar magnetic field B is initially dipolar 
and aligned with the star's spin axis. Then 

i? 3 
B = B p (t = 0) = B d — (2 cos 9e? + sin 9e § ) , (8) 

where B d is the strength of the equatorial magnetic field at the stellar surface and which is observed to be in the range 10 11 — 10 13 
G. The use of a dipolar magnetic field allows us to consider the star as being initially in magnetohydro static equilibrium and 
therefore avoids the calculation of a stable initial configuration. Moreover, we will assume that the electrical currents are 
concentrated at the origin and to avoid singularities we restrict our considerations to a region of the star with radius pR < r < R, 
where < p < 1, and R is the stellar radius. Because any misalignment between the magnetic field and the rotation axis will 
only introduce geometrical corrections 0(1), we expect that all of the features of the generation and evolution discussed below 
will not change significantly when a more generic initial configuration is considered. 

Except for the very first moments after the star's birth, the electrons in the outer layers of a neutron star are strongly degenerate 
and form an almost ideal Fermi-gas. Atoms are partially or fully ionized and form either a strongly coupled Coulomb liquid 
or a Coulomb crystal [Era]. The electrical and thermal transport properties of this dense matter are mainly determined by the 
transport properties of the electrons, which are the most important carriers of electrical charge and heat. At temperatures above 
the crystallization temperature of the ions, the electrical and thermal conductivities are governed by electron scattering off ions. 
Increasingly accurate calculations of the electrical and thermal conductivity of the matter in hot neutron star envelopes can be 
found in the literature ||2q-|2q] and an approximate expression for the electrical conductivity is given by p5]]p| 
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3 Note that expression (H) is roughly correct for densities in the range 10 10 — 10 14 g cm 3 and temperatures in the range 10 6 — 10 8 K, but 
provides a reasonable estimate also at temperatures of ~ 10 9 — 10 10 K which are the relevant ones for the r-mode instability (cf. []28|). 



where T and p are the stellar temperature and mass density. Even for a magnetic field that varies on a length-scale as small as 
AR ~ 0.1R, the magnetic diffusion timescale is 

47r(Ai?)V E „ in6 / AR \ 2 /10 9 K\V P A 3/4 

Tdiff= V " 3X1 ° (lO^j (^-j Uou g cm-3 J y1 '- (10) 

This is well over six orders of magnitude larger than the one year timescale usually discussed for the existence of of r modes 
for a typical newly-born neutron star. As a result, we can neglect the effects of Ohmic dissipation, treating the fluid as perfectly 
conducting on the timescales of interest here. 

B. Lagrangian approach 

In order to simplify our notation, hereafter we will drop the index "1" for the linear velocity so that 6v = 6v±. Combining 
then the equation of mass conservation 

| = -p(V-v), (11) 



where d/dt = (d/dt + v • V), with the induction equation 

OB 



— = V x (v x B) , (12) 



one obtains 
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where we have decomposed the velocity v into v = vo + Sv, with vo = fl x r being the uniform stellar rotation velocity and <5v 
the r-mode velocity perturbation. Here, D/Dt = (d/dt + v • V — fix) is the Lagrangian derivative for a fluid element moving 
at velocity v as seen in the corotating frame. Equation ( |l3| ) can be integrated analytically to give [ ^5p^ ] (see Appendix [§] for a 
derivation) 

— (x,t) = — (x,i )— -^. (14) 

p p ox k (t ) 

Note that V • vo = at first order in the stellar angular velocity and V • <5v = by definition of axial perturbations [cf. 
equations (Q)]. To lowest order in fl and a the flow is therefore incompressible and we can set p(x, t) — p(x, to). The integral 
form ( |l4| ) of the induction equation is particularly advantageous as it shows that the magnetic field at time t and position x can be 
computed from the magnetic field at time to and position x, using the tensorial coordinate strain dx^ (t) / dx k (to) that develops 
between to and t. The advection of magnetic field lines is an obvious consequence of Eq. (pi) and the problem of the magnetic 
field evolution is therefore transformed into the problem of determining the time evolution of the strain tensor Sjk = dx^ /dx . 
While very compact and relatively simpler to solve numerically, Eq. (pi) has the disadvantage of being sensitive to the accurate 
calculation of the strain tensor, which might become difficult when the instability is fully developed. For this reason, and to 
verify the validity of the Lagrangian approach for very large saturation amplitudes, we have also implemented a more traditional 
Eulerian method for the solution of the induction equation which is discussed in the following Section. 

C. Orbit Average Eulerian approach 

We start by rewriting the induction equation ( |12| ) as 

— = V x ((Jv x B) = - ((5v V) B + (B • V) <Jv , (15) 

where now UB/Vt = (d/dt + vo • V — fi x )B = (D/Dt — <5v • V)B is the time derivative of the magnetic field in a coordinate 
system instantaneously corotating with the star. Writing the directional derivatives in (Ob explicitly within an orthonormal 
coordinate system (ef , eg,ei) gives 
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The induction equations (|l6|) are not yet in a useful form since they refer to instantaneous values of the velocity perturbations. 
The magnetic field produced over a single oscillation is simply estimated to be 5B rs -kolB and is uninterestingly small unless 
the seed magnetic field is already very largeQ We need therefore to perform an orbit average of the induction equations ( |16| ) 
over a timescale u~ l <C r <C T, where T is the timescale for a global change of the velocity and magnetic fields. In this case 
we can introduce a time-average (secular) drift velocity v = (6v), whose components at a time t are defined as 

v l (r,9,<p,t) = — I 5v l (r,9,<f>,t + t')dt' , i = 2, 3. (17) 

2r J-r 

and where, of course, v r = 0. For simplicity and because of the smallness of <5B, we can neglect the nonlinear term 
V x ((<5v x <5B)) introduced by the time average and correlating the velocity perturbation with the magnetic field perturba- 
tion. The orbit average of equations ( jig ) then has the only effect of replacing 8v l with v\ 

At this point a number of considerations can be made to simplify the solution of the orbit average of the Eulerian equations 
(|lq). Firstly, we note that the secular velocities must have the same variable separability as the originating instantaneous velocity 
perturbations. As a result, they can be similarly decomposed as [cf. Eq. (fil)] 

v i = f(r)g i (e^)=r e g i (9,cf > ), i = 2, 3 , (18) 

so that dv l /dr = £v l /r, and 8(v l /r)/dr = (I — l)v l /r 2 . By exploiting the property (p"8|) we are therefore able to remove all of 



the radial derivatives in equations ( 16). 



4 Note that although the newly generated magnetic field is small, the magnetic tension forces due to the non-axisymmetric deformation of Bo 
might well be comparable with the driving radiation-reaction force, significantly distorting the character of the r-mode oscillations. See the 
discussion in Sect. IV of paper II. 
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FIG. 2. This panel shows the trajectories 8(t), and <fi(t) of fiducial fluid elements during five oscillations for a mode with £ = 2, a = 0.1, 
and fl = Q. B . Continuous lines refer to fluid elements initially positioned on the (f>o = meridian, while the dotted ones refer to fluid elements 
initially positioned on the <f>o = 7r/2 meridian and have been rescaled as <j>(t) —> 4>(t) — 7r/2 in order to be superimposed on the same plot. 



Secondly, while v^ has a polar dependence, it must be axisymmetric if it is continuous and periodic. This can most easily be 
seen by considering the equations of motion (g). In the case 1 = 1, for instance, </> will be periodic with period n and will change 
sign at <f>o = 7r/2, A similar 7r-periodicity must be expected also for v^ which, however, cannot change sign at <f>o = ir/2, 
nor vary as a function (j>, since either of these two features would make fluid elements converge secularly and thus produce 
discontinuities. This property of v^ is synthesized in Fig. g, which shows the trajectories of fiducial fluid elements during five 
oscillations and the overall displacement that results from these motions (this is the same as the left panel of Fig. pi but here we 
have considered only fluid motions near to the star's equator). In particular, we show fluid trajectories for two different initial 
longitudes, i.e. 4>o = and <f> = it/2, and rescale the latter trajectories so that they can be superimposed on the same plot. As 

expected, <j> has different signs for the two initial data, but in both cases v^ has the same sign, magnitude and polar dependence, 
so that we can set 



dv 4. 



= 



(19) 



Thirdly, while v 9 does not have symmetries, it has some important properties: i) it has an overall periodic (in time) behaviour, 
with period ^> 2-k/ui; ii) the azimuthal dependence is of the form e" n ^; Hi) it is always much smaller than v^; iv) the polar 
dependence involves only the lowest order Legendre polynomials. 
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FIG. 3. Time evolution of the polar profile of v 6 and v^ at the surface of the star and for an i = 2 mode, saturated at an amplitude 
a S at = 10 -2 . Different curves refer to different times and show the velocity increase produced by the growth of the mode amplitude. The 
velocity profiles shown here refer to a time t < 340 s. The velocity and subsequently decreases as a result of the star's spin down (not shown 
here). 



This periodic behaviour of v e guarantees that the polar deformations of the magnetic field are confined to very small scales 
and, when averaged over the relevant timescales, will produce negligible net effects. 

Properties Hi) and iv) are shown in Fig. || where we have plotted the evolution of the polar profiles for v 9 and v^ at the surface 
of the star during the growth of an £ = 2 mode. Different curves refer to different times and show the progressive increase 
of v and v® as a result of the mode's growth. The largest velocity values are reached at saturation (which was here set to be 
«sat = 10~ 2 and occurs at t ~ 340 s) and it should be noted that v^ is almost two orders of magnitude larger than v 9 . Since 
both v 9 and v^ depend linearly on fl [cf. equations |pb], their values progressively decrease after saturation as a result of the 
star's spin-down (This is not shown in Fig. ||.). 

The azimuthal dependence of v 9 is the same that would be imprinted on the magnetic field produced by the poloidal velocities. 
In fact, even an initially axisymmetric magnetic field would acquire non-axisymmetric features driven by v 9 . However, because 
these departures away from axisymmetry are always linearly dependent on v e and on its (/(-derivative, we will assume them to 
be negligible and set 
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As a result of (pH), ( p0| ) and after some regrouping, the induction equations ( |l6J ) take the form 
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Vt r 89 ' 

VB § B f v e v 9 dB 9 B 9 dv 9 



Vf 



iW 



88 



(20) 

(21a) 
(21b) 
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Vt 



00 



Z—i^ r + lv»B' ? -B a v 



Ox* 



COttf 



(21c) 



The set of equations (EJ) has now only time and polar derivatives, and can be more easily solved numerically. Results of 
the numerical integration of (J2J]) and of the comparison with results obtained from the Lagrangian approach will be presented 
in paper II. It is important to underline that the only approximation made in obtaining equations ( pl| ) from the more general 
equations ( jig ) comes from ignoring the deviations away from axisymmetry in the velocity field. In this respect, dynamo theory 
and in particular Braginsky's dynamo model [30|], suggests that such an assumption can only lead to an underestimation of the 
actual growth rate of the magnetic field [ pl| , |l5 ] . 

The Eulerian or orbit average approach presented above is clearly more complicated to implement numerically than the 
corresponding Lagrangian method as it involves the use of a numerical grid on which the set of coupled partial differential 
equations needs to be solved. However, the Eulerian method has also been shown to be more accurate when the 7--mode 
instability has reached saturation and the velocity perturbations at the star's surface are 0(ilR). In paper II we will present a 
close comparison of the two approaches discussed above. 



D. Simplified analytical model of the r-mode instability 



While a detailed discussion of the full numerical calculations will be presented in paper II, in what follows we show how, with 
simple back of the envelope calculations, it is possible to predict the generation of very large toroidal magnetic fields as a result 
of the azimuthal secular velocity field produced by the onset and saturation of the r-mode instability. The estimates discussed 
below will then be confirmed by the full numerical calculations. 

We start by considering the phenomenological expressions for the time evolution of the stellar angular velocity and the mode 
amplitude during the period of activity of the instability [O]. These can be summarized analytically as (see also Fig. 1 of 
paper II) 



a(t) 



r» 

r!o[i + d< +2 (i-i sat )] 

"o exp(t/|r ow |) 

a sa t 



-(2M-2) 



if t < isat i 
if t > teat , 

if t < tsat , 
if t > tsat ■ 



(22) 

(23) 



Here |r GW | is the timescale for the onset of the instability, while fi is me star's initial angular velocity. C\ is a short-hand for 



C x 



647rG(2l + 2)(l-l) 2 ' (1+2 

"^+3" [(2£+l)!!] 2 



21+2 QaL '•" 

1-Qoj 



2 

sat -^0 



f) r 



21+2 



dr 



(24) 



with Q being a nondimensional quantity accounting for the internal structure of the star. For an I = rn = 2 mode [ 12 1 



Q 



pr G dr 



IGirR 2 
Using (|2|) and (|3j), Eq. (g) can be integrated to give 

exp(2t/|r GW |)-l 




pr dr 



(25) 
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if t > t s 



(26) 



Using Eq. (1141), we can now estimate the average magnetic field produced at a given time at the surface of the star as a result of 
the secular drift (fZq). Assuming the initial magnetic field to be predominantly poloidal, i.e. B^(x, to) <S B e (x, to), the toroidal 
magnetic field produced is then 
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B*(x,t) - i? r (x,t )^ + B e (x,i )^^ . (27) 

For a "fiducial" neutron star, with mass M = 1.4M Q , radius i? = 12.5 km, and initial angular velocity il equal to il B , the 
timescale for the onset of an unstable I — m = 2 mode is then |t gw | ~ 37 s. If the initial perturbation has amplitude «o = 10~ 6 , 
the mode will saturate at a sa t = 0.1 after a time i sat ~ 430 s [cf. Eq. (E3h] and the volume averaged toroidal magnetic fields at 
saturation and after one year can be estimated to be respectively 



(AB*(t = i sat )) 



Bq + —B 



n 



i* d 3 x - 3.9 x 10 2 (B%) G , 
(AB*(i = lyr)) ~ 1.2 x 10 8 (B%) G , (28) 



where ((Bq) 2 ) = ((Bq) 2 + (Bq) 2 ) is the initial poloidal magnetic field averaged over the stellar volume V*. Expressions 
show that, in the absence of a back-reaction on the kinematics of the r-mode instability, the toroidal magnetic field is tightly 
wrapped around the star so as to have become about two orders of magnitude larger that the seed poloidal magnetic field in the 
short time necessary for the instability to reach saturation. Moreover, the total magnetic field could be amplified by eight orders 
of magnitude on a timescale of one year. 

More detailed computations of the evolution of the magnetic field will be presented in paper II. However, the simple estimates 
outlined above already show that an initial magnetic field Bq could produce, after one year, an equipartition toroidal magnetic 
field of ~ 10 15 G, i.e. a toroidal magnetic field whose energy is comparable with the rotational energy of the star. Thus, an initial 
magnetic field exceeding 10 8 G, (much below the measured magnetic field in young pulsars) would cause significant departures 
from the standard evolution of the r-mode instability. The impact of these very intense magnetic fields on the existence or 
growth of the r-mode oscillations will be discussed in detail in paper II. There, it will be possible to calculate the strength of 
the magnetic field necessary to significantly distort the first r-mode oscillation, or suppress the instability when this is free to 
develop. 

V. CONCLUSIONS 

We have investigated the onset and growth of the r-mode instability in rotating magnetized neutron stars. Because of the 
high conductivity of the hot neutron star matter and the peculiar nature of the instability which is powered by large scale mass 
currents, it is not possible to ignore the presence of the strong magnetic fields that are expected to accompany newly born neutron 
stars. 

Expanding the perturbed velocity field in powers of the mode amplitude we have derived a second-order analytic expression 
for a secular velocity field which we expect to emerge during the nonlinear growth of the instability. These secular motions 
produce a differential rotation both in the radial and in the polar directions. On an isobaric surface, the secular effects manifest 
themselves as a toroidal drift and couple with any pre-existing magnetic field to produce toroidal magnetic fields that rapidly 
grow in magnitude. In order to study how these kinematic effects interact with a magnetic field, we have discussed two different 
approaches to the solution of the induction equation and have derived sets of equations for the two cases. While we have left the 
discussion of the numerical results obtained to the companion paper II, we have here provided first estimates of the magnitudes 
of the magnetic fields that would be produced as a result of the shearing of a pre-existing poloidal magnetic field into a toroidal 
one. In particular, we have shown that it could be relatively simple to obtain, on the timescale usually discussed for the existence 
of the r-mode instability, magnetic fields that would be in equipartition the rotational kinetic energy of the star. The magnetic 
fields that are produced in this way will influence the r-mode instability either by preventing its onset (when sufficiently strong) 
or by suppressing its saturated development. Precise estimates of the critical magnetic field for prevention and damping of the 
instability will be presented in paper II. 
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APPENDIX A: EQUATIONS OF MOTION: NONLINEAR EFFECTS FROM LINEARIZED EQUATIONS 



In this Appendix we outline the perturbative analysis which allows the periodic and secular parts in the linearized equations 
O) to be distinguished. Here we will assume t = 2, but the formalism presented can easily be generalized to the case of an 
arbitrary L We start by rewriting (g) as 



where 



9(t, 9, 0) = aA sin 9 cos (20 + u>t) , 



4>(t, 9, 0) = -aA cos 9 sin (20 + tot) 



. r\ ( r \ 1 /5! 

A = ui { — Co = u> — —\ — 
\rJ \rJ 8 V tt 



(Ala) 



(Alb) 



(A2) 



Next, we expand the solution of ( Al ) in a series of powers of the mode amplitude a, i.e. we look for solutions of the form 

6(t) = O + aMt) + a 2 Mt) + 0(a 3 ) , 



<f>(t) =0 O + a6 x (t) + a 2 9 2 (t) + 0{a 3 ) . 



Substituting (A3) into dAlh yields 



O(t,0,4>) = aA< [sin0 o cos(0i + 9 2 ) + sin(0i +0 2 )cos9 o ] x 



[cos(20 o + ut) cos(20i + 20 2 ) - sin(20 o + uit) sin(20i + 20 2 )] \ , 



0(i, 0, 0) = -aAi [cos O cos(9 1 + 2 ) - sin(6'i + 2 ) sin# ] x 



[sin(20 o + ut) cos(20i + 20 2 ) + cos(20 o + cot) sin(20i + 20 2 )] 
where we have used the abbreviated notation of 0i for a l 0i and 0; for a l 0i- Making use of the relations 

cos(0 1 + 2 ) = 1 - ^- + 0(a 3 ) , 

sin{0 1 +0 2 )=0 l +0 2 + O{a 3 ), 
cos(20! + 20 2 ) = 1 - 20? + 0(a 3 ) , 
sin(20! + 20 2 ) = 2(0i + 2 ) + 0{a 3 ) , 



equations ( A4 ) take the form 



6(t, 9, 4>)=0o + 9i+9 2 + 6 3 + 0{a 4 ) , 
0(t, 9, 0) = O + 0i + 2 + 03 + 0{a 4 ) 



where 



%{t) = 



l (t) = aA sin 9q cos(20o + cot ) , 



(A3a) 
(A3b) 



(A4a) 



9 2 (t) = a A [9\ cos 9 cos(20 o + cut) — 20i sin 9$ sin(20 o + ujt)] , 



(A4b) 

(A5) 
(A6) 
(A7) 
(A8) 

(A9a) 
(A9b) 

(AlOa) 
(AlOb) 
(A 10c) 
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and 



Mt) = o , 



4>i(t) = —aA cos 8osm(2<fio +u)t) , 



<j>2(t) = oiA {9\ sin O sin(20 o + ut) — 2<j>x cos6>o cos(20 o + uit)] 



(Alia) 
(Allb) 
(Allc) 



Note that, in principle, both a and u> (and therefore A) have a dependence on time which at first order in a produ ces or bits with 
increasing radius of gyration. Considering both a and uj as constant, we can easily find the integrals of equations (AKJ) through 
an iterative procedure. In particular, we substitute into the right-hand-side of the equations for #,(£) and 0j(£), the integrals of 
the equations for 9\_\(t) and 0j_i(t), By doing so we then obtain the following equations of motions 



% = const. 



aA 



9\ (t) = sin 6q sin(20 o + ut) , 



>(*) 



2u 2 



■ sin #o cos 9q sin (20o + ut) , 



and 



const. 



aA 



4>i (t) = cos O cos(20 o + ut) , 



Mt) 



a 2 A 2 

2u 2 



[(2(/>o + ut) (2 cos 2 9 Q - sin 2 6 ) + sin(20 o + ut) cos(20 o + ut) (2 cos 2 9 Q + sin 2 O )] 



(A12a) 
(A12b) 
(A12c) 

(A13a) 
(A13b) 
(A 13c) 



4> 3 (t)=aA 



2<f)\ cos #o sin(20o + ut) — 202 cos^o cos(20o + ut) 



2#i0i sin 6>o cos(20 o + ut) + 8 2 sin 9 sin(20 o + ut) 



(A13d) 



A rapid look at equations ( A.12 ) and ( A13 ) shows that while 9\{t), 9%{t), and 0i(i) are periodic functions, 9%(t) and 2 (t) 
have terms that grow linearly in time and are responsible for an 0(a 3 ) drift in latitude and an 0(a 2 ) drift in longitude. 

As discussed in the main text, using equations (ph to compute the displacement of an element of fluid by expanding 9 and in 
powers of a is not equivalent to considering nonlinear effects in the fluid equations. The key issue is whether the fluid drift given 
by (Oh is at least qualitatively correct. In principle it might not be, because the velocity field obtained by solving the nonlinear 
fluid equations that describe r waves might have additional parts that contribute to the drift and largely or completely cancel the 
drift given by the linear velocity field. Nonlinear solutions for the r waves are not yet available, so it is natural to ask how the 
motions of fluid elements that we have found for the r waves compare with the motions of fluid elements in other, more familiar 
waves. 

Sound waves cause fluid elements to drift as well as to oscillate, as discussed by Landau and Lifshitz pQ], who show how 
to compute the mass current density produced by a sound wave in the Eulerian frame, using the solution of the linearized fluid 
equations. The drift is of second-order in the wave amplitude. Calculating the fluid drift produced by sound waves as we have 
computed the drift produced by r waves, i.e., by using the velocity field (g) in the Lagrangian frame to compute the motions of 
individual fluid elements, gives the drift quoted by Landau and Lifshitz. Shallow water waves are another interesting example. 
Here we outline how the fluid equations can be solved exactly for such waves, using the method of characteristics, and describe 
the exact solution. Consider an infinite train of shallow water waves with wavelength 2irk~ 1 much larger than the depth h of 
the unperturbed water. The evolution of surface waves is then determined by the momentum and mass conservation equations 



d% 



du 



Oh 



dt dx dx 



and -TT- 



dh 

dt 



d(hu) 
dx 



0, 



(A14) 
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where u is the horizontal velocity and g is the gravitational acceleration. Using c = \/gh instead of h, implicit solutions of 
arbitrary amplitude can be expressed in terms of the Riemann invariants, I±, of the rightward-moving (+) characteristics and 
the leftward-moving (— ) characteristics: (dx/dt)± — u ± c. Setting up at t — a right-propagating sinusoidal wave of small 
amplitude y = ac cos(kx), u — 2ac cos(kx), where y = c — c and c = y/gh is the speed of small-amplitude waves, we 
obtain for the right-moving characteristics x(x ,t) — x + c [l + 3acos(kx )]t. As the wave propagates, the leading slopes of 
the crests steepen and form singularities at the earliest crossing of the characteristics, which occurs at t s = l/3ac k. At early 
times (t -c t s ), we can use the above characteristic trajectory to obtain the velocity field to O (a 2 ) from the exact solution. 
The result is u = 2ac cos k(x — c t) + 3a 2 tkc 2 sin2fc(a; — c t). Averaged over the period T = 2tt /kc ot the first term gives 
the drift velocity v^ = 2a 2 c ; the second term describes oscillations at twice the wave frequency, with the linearly growing 
amplitude (3/2)vac kt. Subtracting constant-amplitude oscillations, the trajectory of a fluid element that is initially at x is then 
x = x + Vdt + (3/4)vat cos 2k(x — c t). Calculating the fluid drift produced by this wave in the same way we have done for 
the r waves, i.e., using the velocity field 8v\ in the Lagrangian frame to compute the motions of individual fluid elements, gives 
the same drift velocity as the exact solution, to second order in a. 

This calculation shows that shallow water waves produce a fluid drift and that the drift obtained by using the linear velocity 
field to compute the motion of an element of fluid is correct to second order in a. Stated differently, when the full fluid equations 
that describe shallow water waves are solved exactly, one finds that the lowest-order nonlinear corrections to the linear velocity 
field, which includes all effects of 0(a 2 ), make no contribution to the drift; the drift is given exactly to 0(a 2 ) by the linear 
velocity field. 

This result is so important that it may be worth restating it in a different way. Suppose that we solve the fluid equations 
for sound waves and shallow water waves exactly to 0(a 2 ) and write the solution for the Eulerian velocity field as u(x, t) = 
ocui + a 2 u 2 . When evaluated at a fixed location, neither U\ nor u 2 give a secular drift, i.e., a displacement that is linearly 
proportional to t. u\ is a pure sinusoid and therefore does not give a secular drift, ui is a sinusoid that grows in amplitude and 
therefore does not give a secular drift either, because if at any time we wait one period longer, the displacement from the starting 
point given by u 2 is in the opposite direction and even larger; the displacement given by u 2 is not proportional to t. If now we 
compute the trajectory of an element of fluid by solving dx/dt = u(x, t), we find that the only term that is proportional to a 2 t 
arises from the aui term in the Eulerian velocity field and from the fact that the velocity field at the position of the fluid element 
varies as it moves. The new part of the Eulerian velocity field that arises from solving the fluid equations to second order in 
a does not contribute to the secular drift to 0(a 2 ); the secular drift is already given exactly to this order by the velocity field 
obtained by solving the linearized fluid equations. 

We expect the drift of fluid elements that we have computed from the r-wave velocity field <|3J> to be qualitatively correct and 
perhaps exact to 0(a 2 ). 

APPENDIX B: MAGNETIC FIELD EVOLUTION EQUATIONS: THE LAGRANGIAN APPROACH 



We here briefly show the derivation of the integral form of the induction equation (|l4J) discussed in Section IV B| . Consider a 
Lagrangian coordinate system x. The differential equation (Ob can then be written in terms of a generic vector field A as 



D 

~Df 



(A) = (A • V)<5v 



Equation (Bl ) can be readily integrated for an infinitesimal time interval St from an initial time to, so that 



(Bl) 



A i (t + 5t) = A i (t ) + 



A '^8xl 



6t = AJ ^£j 



a; < (to)+w i (to)<Jt 



^A^t )^-[x\t + St)} 



(B2) 



thereby relating the value of A at the new time t + St, to the value of A at to- The extension of (B2 ) to a finite time interval is 
straightforward and yields 



where Ax^t) = x j {t) - X j (t ) = x 3 (t) - x 3 (t Q ). 



k dx k (t Q ) 



A k (x,t ) 



dx j (t) 

dx k (t ) 



(B3) 
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